home *** CD-ROM | disk | FTP | other *** search
/ Magnum One / Magnum One (Mid-American Digital) (Disc Manufacturing).iso / d18 / nrpas13.arc / SVBKSB.DEM < prev    next >
Text File  |  1991-05-01  |  2KB  |  91 lines
  1. PROGRAM d2r8(input,output,dfile);
  2. (* driver for routine SVBKSB *)
  3. LABEL 10,99;
  4. CONST
  5.    np=20;
  6.    mp=20;
  7. TYPE
  8.    glnparray = ARRAY [1..np] OF real;
  9.    glmparray = ARRAY [1..mp] OF real;
  10.    glnpbynp = ARRAY [1..np,1..np] OF real;
  11.    glmpbynp = ARRAY [1..mp,1..np] OF real;
  12. VAR
  13.    j,k,l,m,n : integer;
  14.    wmax,wmin : real;
  15.    a,b,u : glmpbynp;
  16.    v : glnpbynp;
  17.    w,x : glnparray;
  18.    c : glmparray;
  19.    dfile : text;
  20.  
  21. (*$I MODFILE.PAS *)
  22. (*$I SVDCMP.PAS *)
  23.  
  24. (*$I SVBKSB.PAS *)
  25.  
  26. BEGIN
  27.    glopen(dfile,'matrx1.dat');
  28. 10:   readln(dfile);
  29.    readln(dfile);
  30.    readln(dfile,n,m);
  31.    readln(dfile);
  32.    FOR k := 1 to n DO BEGIN
  33.       FOR l := 1 to n-1 DO read(dfile,a[k,l]);
  34.       readln(dfile,a[k,n])
  35.    END;
  36.    readln(dfile);
  37.    FOR l := 1 to m DO BEGIN
  38.       FOR k := 1 to n-1 DO read(dfile,b[k,l]);
  39.       readln(dfile,b[n,l])
  40.    END;
  41. (* copy a into u *)
  42.    FOR k := 1 to n DO BEGIN
  43.       FOR l := 1 to n DO BEGIN
  44.          u[k,l] := a[k,l]
  45.       END
  46.    END;
  47. (* decompose matrix a *)
  48.    svdcmp(u,n,n,np,np,w,v);
  49. (* find maximum singular value *)
  50.    wmax := 0.0;
  51.    FOR k := 1 to n DO BEGIN
  52.       IF  (w[k] > wmax) THEN  wmax := w[k]
  53.    END;
  54. (* define "small" *)
  55.    wmin := wmax*(1.0e-6);
  56. (* zero the "small" singular values *)
  57.    FOR k := 1 to n DO BEGIN
  58.       IF  (w[k] < wmin) THEN  w[k] := 0.0
  59.    END;
  60. (* backsubstitute for each right-hand side vector *)
  61.    FOR l := 1 to m DO BEGIN
  62.       writeln;
  63.       writeln('Vector number ',l:2);
  64.       FOR k := 1 to n DO BEGIN
  65.          c[k] := b[k,l]
  66.       END;
  67.       svbksb(u,w,v,n,n,np,np,c,x);
  68.       writeln ('    solution vector is:');
  69.       FOR k := 1 to n-1 DO write(x[k]:12:6);
  70.       writeln(x[n]:12:6);
  71.       writeln ('    original right-hand side vector:');
  72.       FOR k := 1 to n-1 DO write(c[k]:12:6);
  73.       writeln(c[n]:12:6);
  74.       writeln ('    result of (matrix)*(sol''n vector):');
  75.       FOR k := 1 to n DO BEGIN
  76.          c[k] := 0.0;
  77.          FOR j := 1 to n DO BEGIN
  78.             c[k] := c[k]+a[k,j]*x[j]
  79.          END
  80.       END;
  81.       FOR k := 1 to n-1 DO write(c[k]:12:6);
  82.       writeln(c[n]:12:6)
  83.    END;
  84.    writeln ('***********************************');
  85.    IF eof(dfile) THEN GOTO 99;
  86.    writeln ('press RETURN for next problem');
  87.    readln;
  88.    GOTO 10;
  89. 99:   close(dfile)
  90. END.
  91.